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Abstract. In this lecture I give a brief review of low-dimensional few- 
body problems recently encountered in attempting a quantitative de- 
scription of ultracold atoms and molecules confined in 2D and ID op- 
tical lattices. Multi-channel nature of these processes has required the 
development of special computational methods and algorithms which I 
discuss here as well as the most interesting results obtained with the 
offered computational technique and future perspectives. 
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1 Introduction 

Recent advances in physics of ultracold quantum gases have opened a unique 
opportunity to study of low-dimensional quantum systems (see, for example 
|ll2j ). In experiments with ultracold atoms and molecules the particle motion 
is "frozen" along one or two variables by using the optical potential and its 
dynamics is limited to a plane or a line. In this case, the physical properties of 
the low-dimensional system can be dramatically altered under the action of the 
optical trap. Moreover, one can effectively manage interatomic interactions by 
varying the parameters of the trap. This opens up new possibilities for studying 
strongly correlated quantum systems [T], few-body and many-body effects in 
low-dimensional systems (2i3) . new mechanisms of molecule formation 4J, and 
possible elements for quantum computer [I]. 

The simplest example of few-body systems, a quantum two-body system, 
acquires in a confined geometry of a ID trap unique properties as compare to the 
conventional two-body problem in a 3D free space [112] . Confining of the motion 
in the transversal directions leads to transverse quantization of the motion and 
makes the pair collisions in a ID trap essentially multichannel ones. Possible 
transverse excitation energies of the particle define the channel thresholds for 
the scattering in the trap. Depending on the initial conditions, the parameters 
of the trap and the interatomic interaction, different multichannel effects may 
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occur at pair atoniic collisions in the ID confining trap. It is a confinement- 
induced resonance (CIR) defined for zero-energy limit when the ground state 
energy threshold of the colliding pair coincides with the resonance in the closed 
transverse channel |5I6) . transverse excitations/deexcitations during the collision 
[7] , center-of-mass effects |8I9I10| , a new mechanism for molecule formation [TT] , 
a CIR-splitting in anisotropic traps [H], anharmonicity effects [TU], and three- 
body effects 

A quantitative description of the mentioned effects has demanded to develop 
new computational methods that we consider in the next section. Some results 
obtained with the elaborated computational techniques are briefly discussed in 
the section 3. Concluding remarks are given in the last section 4. 



2 Methods 

2.1 Quantum Two-Body Problem in Confined Geometry of ID 
Harmonic Trap 

The collisional dynamics of two distinguishable atoms with coordinates ri , 
and masses mi , m2 moving in the harmonic waveguide with the transverse 
potential 5 w^wf {pi = riSinOi) is described by the 4D time-dependent 
Schrodinger equation [h = 1) [918] 



d 

i-Q^^{PR, r, t) = H{pR, v)%Ij{pr, r, t) 

with the Hamiltonian 

H{pR,v) = Hcm{pr) + Hrei{v) + W{pR,v) . (1) 



Here 



and 



z/i ar^ 2 mi m2 

describe the center-of-mass (CM) and relative atomic motions. The potential 
V{r) describes the atom-atom interaction, pR and r = ri — r2 (r, 9, (f)) 
{p, 4>, z) are the polar radial CM and the relative coordinates and M = nii +m2, 
^ = mim2/M. The term^^^^2g# = ___^_^(^ sin 61^ + represents the 

angular part of the kinetic energy operator of the relative atomic motion. 
The term 

W{pr,y) — p{iol — uj2)r PR sin 9 cos (j) (4) 

in the Hamiltonian (1) leads for two distinguishable atoms that feel different 
confining frequencies wi 7^ a;2 to a coupling of the CM and relative atomic 
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motion, i.e. to the nonseparability of the quantum two-body problem in confined 
geometry of the ID harmonic trap. The problem is to integrate the Schrodinger 
equation from time t = to the asymptotic region t — !■ +00 with the initial 
wave-packet 



^{pB.,r,t = 0) = Nr^exp{~^ nJ + ^^o^ (5) 



_P2_ „ - zq ? 

representing two different noninteracting atoms in the transversal ground state 
of the waveguide with at — and the overall normalization constant 

N defined by < ■ip{0)\tp{0) >= 1. We choose zq — >■ — 00 to be far from the origin 
z = and — > 00 to obtain a narrow width in momentum and energy space 
for the initial wave-packet. The wave-packet moves with a positive interatomic 
velocity vq = kof/j = ^j2t^J\i (defined by the longitudinal colliding energy e||) 
to the scattering region at z = and splits up after the scattering into two parts 
moving in opposite directions z — >■ ±cx). If the initial conditions permit open- 
ing inelastic channels it can lead to the collisional excitation of the transverse 
vibrations [7] or formation of molecular bound states [TT]. Multichannel charac- 
ter of the scattering is also developing at the region of the CIR where resonance 
in the closed channel leads to the resonant behaviour of the elastic scattering 
amplitude if the initial colliding energy coincides with the resonant energy of 
molecular state in the closed channel [61917] . 

The modeling of the interatomic interaction Vir) was discussed in [619111] . 



2.2 Discretization of the Angular-Subspace: 2D Nondirect Product 
Discrete Variable Representation 

For discretization of the angular part in (1) we apply the 2D nondirect discrete 
variable representation [13114115] . For ID quantum problems discrete variable 
representation (DVR) was suggested in [16| and broadly applied in the quan- 
tum chemistry computations J[7j. However, it happened that the attempts to 
extend this approach for the case of more higher dimensions as a direct prod- 
uct of ID DVRs were not so efficient [17118 '. In 1997 we have suggested a 2D 
nondirect product DVR (NDDVR) in application to the time-dependent 3D 
Schrodinger equation describing hydrogen atom in a laser field of arbitrary po- 
larization [13'14'. The method happened to be very efficient in application to 
different few-dimensional quantum problems [15I19I20I21] . Mathematical aspects 
of the 2D NDDVR were discussed recently in [18] . 

In the 2D NDDVR we seek for the time-dependent solution r, t) of the 

problem (1-5) according to the expansion [13.15i9j 

AT 

with respect to the two-dimensional basis 

N 

/,(r2) = ^y,(r2)(F-i),, (7) 
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associated with a mesh Slj = {9jg,(f>j^)- For the variable, the Ng mesh points 
9jg correspond to the zeros of the Legendre polynomial Pjve (cos Oj^ ) — 0. For the 
(j) variable, the mesh points are chosen as (j>j^ = 7r(2j0 — 1)/N^. The total 
number N = Ng x of grid points Qj = {Bj^ , (fij^ ) is equal to the number of 
basis functions in the expansion (6) and the number of terms in the definition 

(7) , where the symbol v represents the twofold index {I, m} and the sum over u 
is equivalent to the double sum 

AT {N^-l)/2 Iml+ATfl-l 

E E • (8) 

u=l m=~{N^-l)/2 l = \m\ 

The I and m indexes show the number of zeros over 9 and variables of the 
polynomials Yi,{fl) which we specify in the next paragraph. Due to the definition 

(8) the values N^f, may be chosen only odd. Ng can take on arbitrary values. The 
coefhcients (i^~^)iyj in the definition (7) are the elements of the N x N matrix 

inverse to the matrix given by the values Yj^ = Y^{Qj) of the polynomials 
Yi,{Q) at the grid points flj. It is clear that fj{Qji) — 5jji at such definition (i.e. 
the basis (7) belongs to the class of Lagrange functions [22]) and the coefficients 
V'j (P-Ri ^1 ^) in the expansion (6) are the values of the searching solution ipifR, r, t) 
at the points of angular grid fij: ipj{pii, r, t) = il){pR, r, Qj, t) . 
The polynomials Y^{Q) in Eq.(7) are chosen as 

I' 

where = Sw holds in general, and thus Yy{f2) coincides with the usual 
spherical harmonic with a few possible exceptions for large values of v such that 
we receive the orthogonality relation 

j Y:{Q)Y,, {Q)dQ « ^i^u.Y^'i = (10) 

J 

for all V and v' < N. Here the N weights Xj are the standard Gauss-Legendre 
weights multiplied by 2t:/N^. For most v and v' the above relation is automati- 
cally satisfied because the basis functions Yi,{Sl) are orthogonal and the Gaussian 
quadrature is exact. For these ly we have C; = Sw in Eq.(9). However, in a few 
cases involving the highest I values, some polynomials P™{d) have to be orthog- 
onalized in the sense of the Gaussian quadrature {C\ ^ 5ui for these specific 

1/2 

values of /). With this choice, the matrix Sji, — Xj Yj^ is orthogonal. 

The basis in the form (7) was initially introduced in 23 for the one-dimensional 
angular space — 9. The fact, that values of the searching wave function 
ipipR, r, t) at the grid points of the 2D angular space tp{pR, r, f2j , t) = ijjj {pR, r, t) 
are utilized in the 2D NDDVR, drastically simplifies the calculations |13ll51^ as 
compared to the usual partial- wave analysis. 
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2.3 Splitting-up Method for 4D Time-Dependent Schrodinger 
Equation 

It is an attractive feature of the 2D NDDVR (6-10) that for the grid representa- 
tion fj{f2) the only nondiagonal part of the Hamiltonian (1) is the angular part 
of the kinetic energy operator (see Eqs. (2,3)) 

1 1 /)2 w 



N 



1/2 

which can be diagonalized by the simple unitary transformation Sju — Yj^ 
|13I14I15] . This property has been exploited for developing an efficient algorithm 
with a computational time scaling proportional to the number N = Ng x of 
unknowns in the system of equations [13] 

d ^ 
«^V'j(Pfl, r,t)^Y^ Hjj, {pR, r)il;y {pn, r, t) . (11) 

f 

For the propagation ■^j{pR,r,t„) -> ■ipj{pR,r,tn+i) in time t„ -)■ tn+i = 
tn + At we have developed a computational scheme [5] based on the component- 
by-component split-operator method suggested by G.I.Marchuk in 1971 [23]. The 
Hamiltonian in (11) permits the splitting into the following three parts 

H,,, {pn, r) = hf) {pR) + hf) (r) + U, {pR, p)5,,, (12) 

where 



2 



2M dpi Apl 2Mpl ^A^V ^ 



1 . p^ u}"^ a;2 

Uj{pR,p) = -{miujl + m2Ujl)pl. + —{— + —)p^ + p{ujI - ujl)ppRCOS(j)j . 

Subsequently we can approximate the time-step ipj{pR,r,tn) — > 'ipj{pR,7\tn+i) 
where tn — > tn+i = tn + At according to 

•ip{tn + At) ~ exp(— — Z\tt/) exp i-iAth'-^^) exp (-iZ\</i(°)) X (13) 
exp{-^AtU)i^{tn) + 0{At^) . 
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The time evolution proceeds as follows. For the first and the last steps accord- 
ing to the relation (13) we write the function -0 and the operators exp{~iAtU /2) 
in our 2D NDDVR (7) on the 2D grid {Qj} = {Oje,4'j^}. Since the potential 
Uj{pR,p) is diagonal in this representation the first and last steps represent 
simple multiplications of the diagonal matrices exp(— |Z\if7j (/9_r., p)). Two inter- 
mediate steps in (13) depending on and /i(i) are treated in the basis Y^, (9) 
where the matrix operators U-^^ {pr) and /j'^-* (r) are diagonal with respect to the 
indices m and I. For that we approximate the exponential operators according 
to 

exp(-iZ\<i) « (1 + ^-AtA)-^{l - ^^tA) + 0{At^) , (14) 

which ensures the desired accuracy of the numerical algorithm (13). Thus, after 
the discretization of r (or pn) with the help of finite-differences the matrix 
A possesses a band structure and we arrive at the following boundary-value 
problems 

i - At i 

(1 + -AtA)i;{t^ + _) = (!_ -AtA)ij{tn) , 

which can be solved rapidly due to the band structure of the matrix A. This 
computational scheme is unconditionally stable [H], preserves unitarity and is 
very efficient, i.e. the computational time is proportional to the total number N 
of grid points over the radial and angular variables [H]. The efficiency of the 
computational procedure is based on the fast transformation with help of the 

1/2 

unitary matrix Sj,^ = Xj Yj^ between the two relevant representations: the 2D 
NDDVR (7) and the F^-representation (9). 



2.4 2D Multichannel Scattering Problem as Boundary- Value 
Problem 

In the case of the pair collisions of two identical atoms in harmonic trap the 
problem reduces to scattering of a single effective particle with the reduced mass 
p, off a scatterer V{r) at the origin, under transverse harmonic confinement with 
frequency uj (h—l) [7] 



Here the energy of the relative two-body motion e — e*^'™-* + is a sum of 



V^r) = eV'(r) . (15) 

z,p Z J 

he relative two-body mo 
the transverse e^^'"^^ and longitudinal collision e|| energies. Due to our definition 
of the confining potential, the transverse excitation energies g^"'™) can take the 
possible values gl^"'™^ = g _ = uj{2n + \m\ -I- 1) > of the discrete spectrum of 
the 2D oscillator ^iioj^p^. 

The problem was in integration of the 3D stationary Schrodinger equation 
(15) with the scattering asymptotic conditions 



0«,m(r) = e'''"'=0„,„(p,<^) + /„„'e'''"'l''l(/)„',„(p,(^) (16) 

n'=0 
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at I z |— )■ -|-oo and finding the unknown matrix elements /„„' of the scattering 
ampUtude, which describe transitions between the initial n and final n' channels 
of scattering. The quantum numbers n of the two-dimensional harmonic oscilla- 
tor denote the transversal excitations/deexcitations n n' of colliding atoms. 
In addition to the quantum number n the asymptotic scattering state is defined 
also by the momentum fc„ of the channel 

The energy e defines the number Ug of open channels in p6p . where the longi- 
tudinal collision energy must be positive e|| = e — g^^^'"^^ > Q . 

It is clear that the scattering amplitude depends also on the index m which, 
however, remains unchanged during the collision due to the axial symmetry 
of the problem. Hereafter we consider only the case m = and the index is 
omitted in the following. Due to the axial symmetry of the interactions in (jlSp 
we separate the 0-variable and reduce the problem to 2D one in the variables 
(r,6') (or (z,p)). 

First, we discretize the 2D Schrodinger equation ([T5|) on a 2D grid of angular 
{9j}fj,i and radial {rj}jLi variables. The angular grid points 9j are defined as 
the zeroes of the Legendre polynomial Pno {cos 9) of the order Ng. Using the 
completeness property of the normalized Legendre polynomials which remains 
valid also on the chosen angular grid 

Ng-l 

J2 Piicose,)Pi{cose,,)^XjX,, = Sjj, , (17) 

1=0 

where Xj are the weights of the Gauss quadrature, we expand the solution of 
equation (ITSI) in the basis fj{9) = J2bJ!a^ Pi{cos9){P~^)ij according to 

i=i 

Here is the inverse of the Ng x Ng matrix P with the matrix elements defined 
as Pji — y/Xj Pi{cos 9 j) . Due to this definition one can use the completeness 
relation (flT)) in order to determine the matrix elements {P~^)ij as {P^^)ij = 
^/Xj Pi {cos 9 j). It is clear from (jlSp that the unknown coefficients Uj{r) in the 
expansion are the values 'ijj{r, 9j) of the two-dimensional wave function ?/'(r, 9) at 
the grid points 9j multiplied by \/XjT. Near the origin r — we have Uj{r) ~ 
r — due to the definition and the demand for the probability distribution 
I ip{r,9j) p to be bounded. Substituting (fTS|) into ([15)) results a system of Ng 
Schrodinger-like coupled equations with respect to the A'g-dimensional unknown 
vector u(r) = {Xy^Uj{r)}^'' 



[H(")(r) -t-2(el- V(r))]u(r) = 0, 



(19) 
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where 

H;?M = - ^ E P.'^(' + 1)(P"')'.' ' (20) 

1=0 

V,y W - F(r,0,)%, = {V-(r) + iu;^^^}^^.^.,, ^rsinOj , (21) 

and I is the unit matrix. We solve the system of equations (IT^ on the quasi- 
uniform radial grid [T3] 

rj=R-:^ , J = h2,...,N (22) 

of iV grid points {vj} defined by mapping rj e [0,i? — > +oo] onto the uniform 
grid Xj e [0, 1] with the equidistant distribution Xj — Xj-i = 1/N . By varying 
N and the parameter 7 > one can choose more adequate distributions of the 
grid points for specific interatomic and confining potentials. 

By mapping the initial variable r in Eg. (1191) onto x we obtain 

[W^°\x) + 2{el - Y{r{x))}]u{r{x)) = , (23) 

with 

h(°)(.) = f{x)5,, 7 - ^ P,,/(/ + 1)(P-1),,., , (24) 



where 



The uniform grid with respect to x gives 6-order accuracy for applying a 7-point 
finite-difference approximation of the derivatives in the equation (p3)) . Thus, 
after the finite-difference approximation the initial 2D Schrodinger equation (|15p 
is reduced to the system of N algebraic matrix equations 

3 3 

^ A^_^u,_, + [Aj + 2{6l - V,}]u, + ^ A^+pU,+p = , = 1, 2, _ 3 

u,+af)uj_i-faf )u,_2+af ^u,_3+4^^u,_4 = j - iV-2,iV-l,iV (26) 

where each coefficient A^ is a A^g x Ng matrix, each aj is a diagonal Ng x Ng 
matrix and each gj is a iVe-dimensional vector. Here the functions u_3, u_2, 
u_i and Uq in the first three equations of the system (for j — 1,2 and 3) are 
eliminated by using the "left-side" boundary conditions: Uo = and u j — Uj 
{j = 1,2,3). The last three equations in this system for j = N,N — 1 and 
N—2 are the "right-side" boundary conditions approximating at the edge points 
rAr_2,rAr_i and rjy — R of the radial grid, the scattering asymptotics (ITB|) for 
the desired wave function u(rj). In order to construct the "right-side" boundary 
conditions atj = iV — 2,iV— 1 and N we used an idea of ref.[3S] i.e. 
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the asymptotic behaviour at the edge points rN^2,'rN~i and = R are 
considered as a system of vector equations with respect to the unknown vector 
/„„' of the scattering amphtude for a fixed n. By ehminating the unknowns 
from this system we implement the "right-side" boundary conditions defined by 
Eqs.dini) at j = iV - 2, iV - 1 and N. 

The reduction of the 2D multi-channel scattering problem to the finite- 
difference boundary value problem (pS)) permits one to apply efficient computa- 
tional methods. Here we use the LC/-decomposition [25] or sweep method [27], 
which is a fast implicit matrix algorithm. The block-diagonal structure of the 
matrix of the coefficients in the system of eg nations with the width of the 
diagonal band equal to 7 x Ng makes this computational scheme an efficient one. 

Solving the problem (j26D for the defined initial vector kn and a fixed n from 
the possible set < n < rie we first calculate the vector function ■tp{kmr,6j). 
Then, by matching the calculated vector tp{kn, R,Oj) with the asymptotic be- 
haviour at r = i?, we calculate the n-th row of the scattering amplitude 
matrix describing all possible transitions n — n' = 0, 1, rig. This proce- 
dure is repeated for the next n from < n < Uo and over upon calculating all 
the elements /„„' of the scattering amplitude. 

3 Some Results 

In this section we briefly discuss some of the most interesting results obtained 
with the developed computational methods. 

3.1 Resonant molecule formation in harmonic waveguides 

In our work [11] it was shown that the quantum dynamics of the coupled center- 
of-mass and relative motion for two different colliding atoms in a harmonic trap 
exhibits a resonant formation process of ultracold molecules. With the devel- 
oped splitting- up method (section 2.3), permitting to treat quantitatively the 
collision in the 4D configuration space (pij, r, 0, we have analyzed in detail 
the quantum dynamics of the process and demonstrated the existence of a novel 
mechanism for the resonant formation of polar molecules. The origin of this effect 
is the confinement-induced mixing (4) of the relative and center-of-mass motions 
in the atomic collision process leading to a coupling of the diatomic continuum to 
center-of-mass excited molecular states in closed transverse channels. The pro- 
cess is illustrated by Fig.l where the probability density distribution averaged 
over the angular variables W{pR, r,t) = J \i^{pR, r, 9, (f>, t)\'^{r^p]i)~^ sin9d9d(j) is 
given as a function of time. For wi ^ uj2 (different atoms "feel" different frequen- 
cies wi and in the trap) we observe that a considerable part of the scattered 
wave-packet (see the probability density distribution W^(pij,r, t) aX, t — Sto) is 
located near origin r — after the collision and corresponds to a molecular 
bound state with excitation of the molecular center-of-mass. It is first excited 
state with respect to the center-of-mass motion what is indicated by the node 
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(a) OJi / UJ2 (^) ^1 — <^2 




Fig. 1. Time evolution of the probability density distribution averaged over the 
angular variables W{pii,r,t) — J \'il}{pji,r,9,(j),t)\'^{r^pR)^^ sm9d9d(t). (a) For 
the waveguide with a;i/2.2 — uj2 — 0.02 (b) for uji — uj2 — 0.02. For t — 8to 
corresponding to a time after the collision the insets show a more detailed view 
of W on the scales < W < 8 x 10~^ corresponding to the continuum part of 
the wave-packet. e|| = 0.004, time is given in units of to = ■7t/uj2- 
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of the density distribution W at the point about pji 3. Note, that during col- 
hsion at i ~ 3to the main part of the wave-packet is temporarily in the ground 
state, decays thereafter rapidly into the continuum but part of it goes into the 
excited molecular state. In contrast to this the case wi = U2 in Fig. lb (the case 
of the center-of-mass separation) shows an almost complete decay the molecular 
ground state into the continuum: The remaining minor part of the probability 
density distribution near r = representing the molecular part is much smaller 
compared to the case uji ^ i02- 

It was also shown that the molecular formation probabilities can be tuned 
by changing the trap frequencies wi and 102 characterizing the transverse modes 
of the atomic species in the trap. 

3.2 Multichannel Scattering and Confinement-Induced Resonances 




Fig. 2. The total transmission coefficients T for bosonic collisions as a function 
of the dimensionless energy s = {e — uj)/uj for the two cases of the system 
being initially in (a) the ground transverse state n = and (b) the first excited 
transverse state n — 1, for several ratios of as/a± and w — 0.002. The black 
solid curve corresponds to as/a± = 0.68 for which the zero-energy CIR in the 
single-mode regime is encountered [S]. 



The grid method for integration of the 2D multichannel scattering problem 
as a boundary- value problem (section 2.4) was applied to extensively analyze 
the transverse excitations and deexcitations as well as resonant scattering pro- 
cesses at pair atomic collisions in harmonic waveguiedes jj . Collisions of identical 
bosonic and fermionic as well as distinguishable atoms in harmonic traps with 
a single frequency uj permitting the center-of-mass separation were explored in 
depth. In the zero-energy limit and single mode regime the well-known CIRs 
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T=0.24 

T=0.67 




Fig. 3. The probability density \tp{x, z)\'^ for bosonic collisions as a function of a; 
and z a.t as/a± = +4.39 for two cases of the single-mode regime (left subfigures 
with e < 1) and two- mode regime (right subfigures with 1 < e < 2). The 
corresponding transmission values T are also indicated. All subfigures are for 
cjO = 0.002 and n = 0. The inset shows a more detailed view of \tjj{x,z)\'^ near 
the origin z = x = 0. 
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for bosonic [5^, fermionic [55] and heteronuclear [8129] collisions were repro- 
duced. In case of the multi-mode regime up to four open transverse channels 
were considered. Important scattering observable in the quasi-lD scattering is 
the transmission coefficient T'„(e) = X^n' ^ I ^nn' + fnn' 

Fig. 2 illustrates the behaviour of the transmission T with varying energy e 
(permitting opening of the four lowest channels during the collision) and scat- 
tering length as/a± (here is the scattering length in the 3D free space and 
a± = y/l/{fj,uj) is the width of the trap). In Fig. 3 the probability density dis- 
tribution I tp{x,z) p calculated as a function of the the transverse (x) and 
longitudinal (z) variables is given. Here the case T = illustrates the appearing 
of the CIR in the scattering. The dual CIR |8I29| leading to a complete quantum 
suppression of atomic scattering was also analyzed in multi-channel regimes and 
possible applications were discussed. 

4 Conclusion 

In this lecture the computational methods have been considered which were 
elaborated for a quantitative analysis of multichannel effects in low-dimensional 
few-body systems. The efficiency of the methods was demonstrated in applica- 
tions for different quantum dynamics of two-body atomic collisions in harmonic 
ID traps. It is very promising further development of these approaches in appli- 
cations to hot problems of the few-body physics in confined geometry of optical 
traps. Thus, the description of the effects of anisotropy and anharmonicity of the 
traps as well as strong anisotropy in the interparticle interactions (scattering of 
polar molecules) requires an extension to the increasing number of spatial vari- 
ables. Progress in this direction also opens up new possibilities for a quantitative 
analysis of more complex low-dimensional few-body systems. 
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